2_Control/6_Discrete_control_and_analysis/Simple discrete controller simulation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simple discrete simulation: Dahlin controller\n",
    "\n",
    "This notebook replicates Figure 17.11 in Seborg et al without using analytic methods. If you want to get some insight into the math behind the z transform, applied to the same problem see the [Dahlin controller](Dahlin%20controller) notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import tbcontrol"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "tbcontrol.expectversion('0.2.1')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from tbcontrol import blocksim"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy\n",
    "import control"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simple discretisation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The control module supplies an easy access point to a transfer function which represents s."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "s = control.TransferFunction.s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "K = 1\n",
    "τ1 = 5\n",
    "τ2 = 3\n",
    "θ = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "Δt = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert θ % Δt == 0  # this is only correct when the delay is a multiple of the sampling time\n",
    "N = int(θ/Δt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "G = K/((τ1*s + 1)*(τ2*s + 1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Discretize G assuming a zero order hold in front of it, this corresponds to using Table 17.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "z = control.TransferFunction.z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\frac{0.02797 z + 0.02341}{z^2 - 1.535 z + 0.5866}\\quad dt = 1$$"
      ],
      "text/plain": [
       "\n",
       "  0.02797 z + 0.02341\n",
       "----------------------\n",
       "z^2 - 1.535 z + 0.5866\n",
       "\n",
       "dt = 1"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Gz = G.sample(Δt)*z**(-N)\n",
    "Gz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Do discrete transfer function math"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dahlin controller relationship"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "λ = Δt   # this is a tuning parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = numpy.exp(-Δt/λ)\n",
    "Gdc = 1/Gz * (1 - A)*z**(-N - 1) / (1 - A*z**-1 - (1 - A)*z**(-N - 1))  # eq 17.65"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After these calculations Gdc is"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\frac{0.6321 z^4 - 0.9705 z^3 + 0.3708 z^2}{0.02797 z^4 - 0.004556 z^3 - 0.02341 z^2}\\quad dt = 1$$"
      ],
      "text/plain": [
       "\n",
       "  0.6321 z^4 - 0.9705 z^3 + 0.3708 z^2\n",
       "----------------------------------------\n",
       "0.02797 z^4 - 0.004556 z^3 - 0.02341 z^2\n",
       "\n",
       "dt = 1"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Gdc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that there are \"extra\" orders of z above and below the line in that expression - we could simplify it. This is what `.minreal` does."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "Gdc = Gdc.minreal()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\frac{22.6 z^2 - 34.7 z + 13.26}{z^2 - 0.1629 z - 0.8371}\\quad dt = 1$$"
      ],
      "text/plain": [
       "\n",
       "22.6 z^2 - 34.7 z + 13.26\n",
       "-------------------------\n",
       " z^2 - 0.1629 z - 0.8371\n",
       "\n",
       "dt = 1"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Gdc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Much better."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Convert from positive to negative powers of z\n",
    "\n",
    "We still have a problem, though, since the control library uses positive powers of $z$ rather than negative powers of $z$. `tbcontrol.conversion` contains some methods to help convert such polynomials. Also notice that the `.num` and `.den` properties of the control library's tf objects are designed for MIMO systems, so they are nested lists corresponding with the index of the output and the input we want. Since we're just looking at one input and one output, we need the first element of the first element (`[0][0]`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "import tbcontrol.conversion"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "Gdc_neg_num, Gdc_neg_den = tbcontrol.conversion.discrete_coeffs_pos_to_neg(Gdc.num[0][0], Gdc.den[0][0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[22.59988127611586, -34.69674036625465, 13.258134912008916]"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Gdc_neg_num"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1.0, -0.1628886995509321, -0.8371113004490679]"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Gdc_neg_den"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simple blocksim simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "blocksim_G = blocksim.LTI('G', 'u', 'yu', G.num[0][0], G.den[0][0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "blocksim_Gdc = blocksim.DiscreteTF('Gc', 'e', 'u', 1, Gdc_neg_num, Gdc_neg_den)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "diagram = blocksim.simple_control_diagram(blocksim_Gdc, blocksim_G, ysp=blocksim.step(starttime=5))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "ts = numpy.arange(0, 20, 0.01)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "result = diagram.simulate(ts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_outputs(result):\n",
    "    fig, (ax_y, ax_u) = plt.subplots(1, 2, figsize=(15, 5))\n",
    "\n",
    "    ax_y.plot(ts, result['ysp'])\n",
    "    ax_y.plot(ts, result['y'])\n",
    "    ax_y.set_ylim(0, 2)\n",
    "    ax_y.set_ylabel('y')\n",
    "\n",
    "    ax_u.plot(ts, result['u'])\n",
    "    ax_u.set_ylim(-50, 50)\n",
    "    ax_u.set_ylabel('u')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA34AAAEzCAYAAABucNGPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd5xdVbn4/8+TmRSSQApJKCkkhNCLQAARlKY0KfYLlos113q5+r3X3i7qtXexcBFFfwooisZrEKOAYAEJSG8JoSQkkp5AQsrMrN8f+0wYhplkZvbec8p83q/XeZ1zdjvP7JycfZ6z1npWpJSQJEmSJDWuQdUOQJIkSZJULhM/SZIkSWpwJn6SJEmS1OBM/CRJkiSpwZn4SZIkSVKDM/GTJEmSpAZXWuIXEZMj4rqIuC8i7omI87vYJiLiGxGxICLujIjDOqw7LyLmV27nlRWnJEmSJDW6KGsev4jYDdgtpXRbROwI3Aq8LKV0b4dtTgfeA5wOHAV8PaV0VESMBeYBM4FU2ffwlNLqUoKVJEmSpAZWWotfSmlpSum2yuMngfuAiZ02Oxv4UcrcBIyuJIynAHNTSqsqyd5c4NSyYpUkSZKkRtYvY/wiYipwKHBzp1UTgUUdni+uLOtuuSRJkiSpl5rLfoGIGAn8AviPlNK6zqu72CVtY3lXx58FzAIYMWLE4fvuu2+OaCVJ9eDWW29dkVIaX+046sW4cePS1KlTqx2GJKkfdHeNLDXxi4jBZEnfT1JKv+xik8XA5A7PJwFLKsuP77T8+q5eI6V0EXARwMyZM9O8efNyxy1Jqm0R8Wi1Y6gnU6dOxeujJA0M3V0jy6zqGcD3gftSSl/pZrPZwL9Wqns+H1ibUloKXAOcHBFjImIMcHJlmSRJkiSpl8ps8TsGeANwV0TcXln2YWAKQErpu8AcsoqeC4ANwJsq61ZFxKeAWyr7XZBSWlVirJIkSZLUsEpL/FJKf6brsXodt0nAu7pZdwlwSQmhSZIkSdKA0i9VPSVJkiRJ1WPiJ0mSJEkNzsRPkiRJkhqciZ8kSZIkNTgTP0mSJElqcCZ+kiRJktTgTPwkSZIkqcGZ+EmSJElSgzPxkyRJkqQGZ+InSZIkSQ3OxE+SJEmSGpyJnyRJkiQ1OBM/SZIkSWpwJn6SJEmS1OBM/CRJkiSpwZn4SZIkSVKDM/GTJEmSpAZn4idJkiRJDc7ET5KkGhYRTRHxj4j4v8rzaRFxc0TMj4grImJItWOUJNU+Ez9Jkmrb+cB9HZ5/HvhqSmkGsBp4S1WikiTVFRM/SZJqVERMAl4KXFx5HsCJwJWVTS4FXlad6CRJ9cTET5Kk2vU14P1AW+X5zsCalFJL5fliYGI1ApMk1RcTP0mSalBEnAEsSynd2nFxF5umbvafFRHzImLe8uXLS4lRklQ/TPwkSapNxwBnRcQjwOVkXTy/BoyOiObKNpOAJV3tnFK6KKU0M6U0c/z48f0RrySphpn4SZJUg1JKH0opTUopTQXOAa5NKb0OuA54VWWz84BfVylESVIdKS3xi4hLImJZRNzdzfr/iojbK7e7I6I1IsZW1j0SEXdV1s0rK0ZJkurQB4D3RcQCsjF/369yPJKkOtC8/U367IfAt4AfdbUypfRF4IsAEXEm8N6U0qoOm5yQUlpRYnySJNWFlNL1wPWVxwuBI6sZjySp/pTW4pdSugFYtd0NM+cCl5UViyRJkiQNZFUf4xcRw4FTgV90WJyA30fErRExqzqRSZIkSVJjKLOrZ0+dCfylUzfPY1JKSyJiAjA3Iu6vtCA+RyUxnAUwZcqU8qOVJEmSpDpT9RY/skplz+rmmVJaUrlfBlzFNsYyWK5akiRJkratqolfRIwCjqNDKeqIGBERO7Y/Bk4GuqwMKkmSJEnavtK6ekbEZcDxwLiIWAx8AhgMkFL6bmWzlwO/Tymt77DrLsBVEdEe309TSr8rK05JkiRJanSlJX4ppXN7sM0PyaZ96LhsIXBIOVFJkiRJ0sBTC2P8JEmSJEklMvGTJEmSpAZn4idJkiRJDc7ET5IkSZIanImfJEmSJDU4Ez9JkiRJanAmfpIkSZLU4Ez8JEmSJKnBmfhJkiRJUoMz8ZMkSZKkBmfiJ0mSJEkNzsRPkiRJkhqciZ8kSZIkNTgTP0mSJElqcCZ+kiRJktTgTPwkSZIkqcGZ+EmSJElSgzPxkyRJkqQGZ+InSZIkSQ3OxE+SJEmSGpyJnyRJkiQ1OBM/SZIkSWpwJn6SJEmS1OBM/CRJkiSpwZn4SZIkSVKDKy3xi4hLImJZRNzdzfrjI2JtRNxeuX28w7pTI+KBiFgQER8sK0ZJkiRJGgjKbPH7IXDqdra5MaX0vMrtAoCIaAIuBE4D9gfOjYj9S4xTkiRJkhpaaYlfSukGYFUfdj0SWJBSWphS2gxcDpxdaHCSJEmSNIBUe4zf0RFxR0RcHREHVJZNBBZ12GZxZZkkSZIkqQ+aq/jatwF7pJSeiojTgV8BM4DoYtvU3UEiYhYwC2DKlCllxClJkiRJda1qLX4ppXUppacqj+cAgyNiHFkL3+QOm04ClmzjOBellGamlGaOHz++1JglSZIkqR5VLfGLiF0jIiqPj6zEshK4BZgREdMiYghwDjC7WnFKkiRJUr0rratnRFwGHA+Mi4jFwCeAwQAppe8CrwLeEREtwNPAOSmlBLRExLuBa4Am4JKU0j1lxSlJkiRJja60xC+ldO521n8L+FY36+YAc8qIS5IkSZIGmmoWd5EkSRrQ1m7Ywu2L1+Q+zsTRO7DXhJEFRCSpUZn4SZIkVcn/zLmPK+Yt2v6G27HTsGbu/OQpBUQkqVGZ+EmSJFXJU5tbmDh6B75x7qF9PsZPbnqU39zZbQF0SQJM/CRJkqpq2OBBHL7HmD7v/8f7nigwGkmNqmrTOUiSJA14CSqzW+U7TCogFkkNzcRPkqQaFBGTI+K6iLgvIu6JiPMry8dGxNyImF+573tTkaoukT9jKyBvlDQAmPhJklSbWoD/l1LaD3g+8K6I2B/4IPDHlNIM4I+V56pTKUEReZsNfpK2x8RPkqQalFJamlK6rfL4SeA+YCJwNnBpZbNLgZdVJ0IVxRY7Sf3BxE+SpBoXEVOBQ4GbgV1SSkshSw6BCdWLTHllLX75Mr+8+0saGEz8JEmqYRExEvgF8B8ppXW92G9WRMyLiHnLly8vL0DlUsQYP4BkdRdJ22HiJ0lSjYqIwWRJ309SSr+sLH4iInarrN8NWNbVvimli1JKM1NKM8ePH98/AatP7OopqT+Y+EmSVIMiq/H/feC+lNJXOqyaDZxXeXwe8Ov+jk3FKaKhzsRRUk84gbskSbXpGOANwF0RcXtl2YeBzwE/i4i3AI8Br65SfCpAUR007egpaXtM/CRJqkEppT/TfaX/k/ozFpUnFTSBuyRtj109JUmSqihv2mfaKKknbPGTJEmqmlTIGL28YwVTSnzu6vtZtHpDruMEwRuPmcoRU8fmC0hS4Uz8JEmSqqSQWRgKyByf2tTC925YyLiRQxgzfEifj/PQ8qcYO2KIiZ9Ug0z8JEmSqqiWhvi9/bjpvPWFe/Z5/8M+NbewuQklFcsxfpIkSVWSyLpHVpupmtT4TPwkSZKqJBXQ17P6aeMzgoK6r0oqnImfJElSlSSK6+pZRBIpqXGZ+EmSJFVRLbTYFZUzRthtVKpVJn6SJElVkgpo8iuyOIyTyUuNy8RPkiSpSopsHcvValdYIOEYP6lGmfhJkiRVSUqpJrp6Smp8pSV+EXFJRCyLiLu7Wf+6iLizcvtrRBzSYd0jEXFXRNweEfPKilGSJKna8vauLGI6iPa59/IeKftbbPKTalGZLX4/BE7dxvqHgeNSSgcDnwIu6rT+hJTS81JKM0uKT5IkqWHUQrpl66VUu5rLOnBK6YaImLqN9X/t8PQmYFJZsUiSJNWilIpqZStGEcdyjJ9Um2pljN9bgKs7PE/A7yPi1oiYVaWYJEmSSldUJc088/gVOZ2DpNpUWotfT0XECWSJ37EdFh+TUloSEROAuRFxf0rphm72nwXMApgyZUrp8UqSJBUl1UQHzWcUkbfZ4ifVpqomfhFxMHAxcFpKaWX78pTSksr9soi4CjgS6DLxSyldRGV84MyZM/2okSRJdaOQrp5FxFHAMaCYQjOtbYn7lq6jtS1fVDvtMJhp40bkjkdqFFVL/CJiCvBL4A0ppQc7LB8BDEopPVl5fDJwQZXClCRJKk1KxXWPrJVfv/O2Yl556yI+8Iu7ConlxvefwOSxwws5llTvSkv8IuIy4HhgXEQsBj4BDAZIKX0X+DiwM/DtSt/2lkoFz12AqyrLmoGfppR+V1ackiRJ1VREK1lR8o43jMjf1XPd0y0AXPjaw9hhSN/KUdz26Bq+dd0C1m3cki8YqYGUWdXz3O2sfyvw1i6WLwQOee4ekiRJjaWIMX7FVOKslfbCZ87JcfuMZ+TQvn1V3dKaHaOG/iyp6mqlqqckSdKAkxKFTX5XRJKTfzL5/F1O2/+OPKHUThuqVDtM/CRJkqqkwLwvl1pqGCsilqKmyJAaiYmfJElSHaulJCciamJOwPZd7eopPcPET5IkqVoKrepZwHjBAuLI65munvmjqbV5EqVqMvGTJEmqkkSqiaqeRbaM5U22iix4Y4uf9AwTP0mSpCoqrMWvmMFx1dy9sGNtTfyKCUVqCCZ+kiRJVVIrLVKFdonMeahCqpNWWlFraZoKqdpKm8dPkiRJ25YoYAqFIlvZ8u5f/V6rmYJa/H7wl4e5cf6K3OGcvP8unHPklNzHkfIw8ZMkSaqSlGpjjF+RimpjK6KqZ14/vulRVjy5iT12HtHnYzyyYj0rntpk4qeqM/GTJEka6IqagqGISpyV7pl5jtU+xUXenp4pwXH7TOCb5x7a52O86Qd/Z8VTm/MFIhXAMX6SJElVUkhXzwJbDIvoqpl3XF2xw/LyxpIK6P4aTiuhmmDiJ0mSVCWFTqOQ41hFds+spa6euVv8csYBMChqp4iPBjYTP0mSpCrJEovaGeNXC+MN23OkPJEUNZ1DSkWMFwzaTPxUA0z8JEmS6lgN5Y0ExYyryx9HQWP8SLkT86zFz8xP1WdxF0mSpGopYAzZ1kPlaN8qKi8psvUyz7G2tvgVMN6wiCku8p7fRas2cNrXb+SpTS25jjOkaRA/fPMRvGD6uHwBqS6Z+EmSJFVJMcVdilNIcZfc+7dX9ey7wpLplP9gQf7iLo+veZqnNrXwikMnMmns8D4dY/2mFr7/54d5dOUGXjA9VziqUyZ+UhkW/R0e/QsMHwf7nQk7jK52RJKkBpevuEsxTX5FJFztf0ctJKGQf9zjoEH5W/zaKgd49czJHD195z4dY9mTG/n+nx/eeiwNPCZ+UpFaNsPsd8OdVzyz7A+fhFf+L0w/sWphSZJqUzHFQ4pTTOJWA4nF1q6e+Q6TUmJQAS1+uZOtyu55YhlUyaTbrDQzYJn4SUVJCX71Drj7SnjR++Hod8HKh2D2e+Cn/wJvuAqmHlvtKCVJNaSI4iG1VNyliMxxa1XPPGP82ou75J3HjwLObwFTXLRtbQXtezBbE7+cwVx911Iu+cvD+Q4C7DVhRz77ioNyH0c9Z1VPqSi3/zRL+k78KJz4kax756TD4U2/hdF7wJVvhqeW9V88C/8El70WvrQ3fGV/+MVbYfGt/ff6kqTtKnQevzz71kgc2QHyB7M1Pyqgwmjurp4FTG7Y3mKYr8Xv2cfqq2vu+Sd3Ll7L4KZBfb79c91GLvv7Y7XROjyA2OInFWHjOvj9R2HKC+DY9z173Q5j4DWXwkXHw+8/Bq/4XrmxtLXC1e+HWy6GkbvCXi+G1s0wfy7c9XN4wXvgxf8Ng5rKjUOS1CM11WBXQ4Vm8igo76MtpULOSd5kq5BW0IJa/FoT7D56B376tuf3+Rjf/ON8vjz3QVrbEs1Nff+bnt7cypMbt/R5fwACxo8cWlPzaZbFxE8qwt8uhKdXwan/03VCtcsBWcJ145fhiLfC5CPKiaO9u+mdV8DR74aTPg7NQ7N1G9fBHz4Bf/0mrH0cXnmxyZ8kVVlKRSRb+b+wFtruUsBwttznJIqax6+IWIpJQKGYFr+8rWxtbfnHPQ6qHKA1pT4nIy2tbRzz+WtZtX5zvmCA80+awXtfsnfu49Q6Ez8pry1Pw83fhX3PgN0P7X67Y98Ht14Kf/o8vP7KcmL527eypO+Ej8Bx73/2umE7wRlfhdFTsoIzIyfAaZ8vJ46OVj2ctTQuvgWeXgMjxsGkI+DAV8CYqeW/viTVsPwTF3Q4VhFdJHPGUkSrSVFz50EBY/wS5I1mUEQhRWYg3/ltGtTe4pcz8Utp67HyxtKao/lxc2sbq9Zv5tQDduWFe/d9XsIv/O4Blqx5us/7A6x9egu/vXMpLW1tuY4D8Nojp9DcVM5oPBM/Ka+7fwkb18BRb9/2dkNHwtHvhD9eAEvvgN0OKTaOJ+6BuZ+A/c6CF/1X99sd+95srOFN34Ypz4cDXl5sHO02PZl1f73tR9mVc8L+MGLnrODNA3Pg2k/DIefCyZ+C4WPLiUGS1K+C4qaGqA210dWzPZ8ooqpna87cpLUtbT1WXzUXkPhtac32PWLaWF531B59Ps73/rSQLTlPyuzbH+djv74n1zHavWbmZJpL6pBl4iflNe8SGLdPzyp2HvFW+PPXsu6Wr7y4uBhSgt/+Pxg2Cs78+vb7pbzkgmyuwdnnw5SjYcddi4sFYPWj8OOXw+qH4ch/g2POh512e2b9mkVw03fg79+DBXPhVT+AqccUG0NnmzfAI3+GJf+Ap54AEoyYALseBHu8wORTUlWkIsaQFVFJc+vkefmPlVchlU7bj1VIcZecopgup0CuhKt91yJa/PImfs8koX2PpX3f5pytj4Obgi05Bz6u39wKwJ8/cAI7DM6XtQ1tLq/2ZqmJX0RcApwBLEspHdjF+gC+DpwObADemFK6rbLuPOCjlU0/nVK6tMxYpT5ZsQAenwen/E/PrrzDRsHB/5K1gp22qrhk477fwGN/gzO/0bNjNg2GV1wE3z4arvkwvOqSYuIAWP0IXHIabNkA5/1f1wnd6MnZeMhDzsmqnf74ZfDy78KBrywujnbrlsKNX4I7roDNTwLxzDnasApIMKgZZpySTcFRdgLabsOqrOV183poHgI7jIUdd8tm+pW2IyJOJbt+NgEXp5Q+V+WQlENRuVYttLNFEUlOoV09c8ZC/sQ6b5IE+ZO1jnHkHuOXyN3Vs72gS57Er71bZZ7iMACDmwaxpSVfi9/myv677jSstG6aRSi7xe+HwLeAH3Wz/jRgRuV2FPAd4KiIGAt8AphJ9n/u1oiYnVJaXXK8Uu/c+6vsfv+X9Xyfw8+DW/4X7vwZPH873UN76m8XZlNGHPr6nu+z8/Ss2+efPgeHvgGmn5A/jtYtWSK3ZT286eqsqM227HYwvHVuNu3EL94GTUNhvzPyx9Huth/DnP+CthY4+DVw0Ktg0pFZt1uALRuzFsAHfgt3XA4//C1MPwlO+wKM26u4OCBL8O77DTz4O3j0b/DUP5+7zeAR2Tnb8zjY+zSYeFi5E3RtXAdrF8OTS7Jz0boZBg/PKtGOHJ+9pywAVHMiogm4EHgJsBi4pXKNvLe6kakvaqWafaHTOdTE31RMkpNSyj/ukQKqem4t7lL9efxaiyjuUkCLX0trMS1+zU1BS86TsrmljUFBTSd90IPELyLeDfykL0lXSumGiJi6jU3OBn6UsnfzTRExOiJ2A44H5qaUVlVimAucClzW2xikUt37qyyRGDWx5/vselBWBOa2S+Gof8v/xX7xrbDoJjjls73/kn7se+HOy+F3H4S3/wWacv4WdO2n4fFb4dWXbj/pa7fDGHjdz7NWvyvfBK//BUx7Ub442toqFUy/AXseD2d8DcZOe+52g4fBHkdntxM+kk2BccMX4bvHwEmfyMZt5m2B27AK/vwVuPVHsGktjNwli2mXA7P3zeAR0LoJ1q+AFQ/C47dl1V9v+CLsvBfMfDM877XZecprw6psfOXC67Ouvmse3fb2TUOzGCYdDnscm7WGjpqUP44GFBEf72p5SumCEl7uSGBBSmlh5bUvJ7uemvjVqaJ+3yki4co/V3kxFUaL6v5aCy1+RbWCQr5LUlHz+LWltLUqZ181dajq2VftSWNTzuv04KZBucf4bW5tY3CNJ33Qsxa/Xcl+TbwNuAS4JhU32+JEYFGH54sry7pbrjp075J1vPeK29mcdzRxjdm9bSk/2XAX3x7yZn7+pet7te+ZW47ifZu+y7994fs82JSvZemjG7/EUQznX/48lQ1/6V0cAMe0nMunN36Wb3z+g1w1pO+tbYe33M6XNn6N3zSfwleuHg1X9y6WHdP5fD19mPGX/gvn7/BZFjZN7VMcw9JGPrzxK7yw9WZ+3Xwa33jibbRd8iiwnSQHgIMYO+hr/L/0bV5wzYe47Q+X8dmh57NiUO+rhQ1OW3jZljm8YfMVDOdp/tT8Aq7a4aXck/YlLRwEC58VNTAKmA6cxk7D13F0yy2csfr3HHjNh3n6mv/m980n8MshZ/DYoMm9imNEWs+xLTdxYsuNHNZ6J820siLGcM+gfXlgyIv4Z+zC8kE78zTDaKWZoWxip/Qk49NKJrctZurKRRy47BeMvC3ruLE4duO2poO5rfkQ7hm0Dyti5+1+K9ohbWB626Ps1bqQGW0LGXbQWZx49nm9+jvqwPoOj4eRDXO4r6TX6uoaeVRJrwXAdfcv44L/M68sw+LVG5g2bkS1wyhUEZU0i2hlqwSTO5YixrPlPSftDVL5xvgFEdl0DHm0tiWacp6T9sSvvdWuL9pb6Qbn7eo5qIDEr6WNISWOzSvKdhO/lNJHI+JjwMnAm4BvRcTPgO+nlB7K+fpd/UulbSx/7gEiZgGzAKZMmZIzHJXhniVreeCJJ3nxfhMYPqRx6gkds+Z62AArJ57AQUNG9WrfJ1pfypaHvs9rh/6FX+xyeJ9jGL1lGccv/At/GvMqpk/Ybfs7dGFdOpkHFl/DWzdextLJL2V98+heH2Nky2o+8MjXWTpkKtfv8V4OGjSsD5GM4gdbvsT7HnsnX97yKb6624WsHrxLr44wessy3vb4R5nY+hBXTvh3bhj9cg7o9cVpFJenL/DI2jm8ctk3uXTj+cweP4ubR51Gawze7t7NbZs59MnrOHXlpYzfsoT7hh/Brya8g6VD9ySA5wx27iaGJUzmIl7BxI3zedGaqzh93VzO3vA7Hhu6N/eOPIpHh+3H0qHTWNe0My2DhkBKNKctjG5ZxoTNi5my8QH2evp29nz6bprTFlYO3pXrRr2G23c8nkVD935WshbA8A6vvpEsq1gE/BWI1MrumxYyY8Pt7L3hNk5++kbO2ngNAJtiGCsG786Gph3Z2DSSFpppTlsYnDazY8sqRresYETbuq3HXj9oJx5NR/foLNSTlNKXOz6PiC8Bs0t6uR5dI4u8Po4aPpiDJvbuc049c9DEUbzy8Hwt6UVOPJ27qEoRhWYKHK2YPwktYIqMyN+9sr2VLu/pHRRRSCy5W/y2djvN09UzS9byjjccOngQN85fwYyPzOl7LG2JcSOH5oqjP/ToW3hKKUXEP4F/Ai3AGODKiJibUnr/tvfepsVAx5+vJwFLKsuP77T8+m5iuwi4CGDmzJk10atcXfvEmQcweezw7W9YhCf/mY1XGl3ijwE//zps2J2P/euZfbvS/eylHPfw9Rz36u9lxT36Yu7H4WE48V8/yolj+l7KmCcuhO8ey2fH/h+89Mvb376jtjb46ashNjDqLf/Hl3vaxbPbWGbDJafy36s/AK/9OYzv4YSqD98IV74L0tPw2it41d4n86pcgRwGK8+FX7+Lcx77Cuds/Hk2/cS+p2fdNJs6JIEb12VjBR/8Hdx1Jaxflm3zkm+y314vZr9ccRwKvAbWr4R//JgpD1zNlMU/gdTh18mmoZBas7GMW0XWrfiQf4P9X8bOk2bykghe0uc4ZmZxALS2ZH/v0tsZumohE1c/AhvXZuehdXP2fm4aCiP3zQrWjJoIEw6AXQ9kxE4T2b/McYu1YziwZ0nH7u7a+SxFXh8PmzKGw6YU0N1Y5aqRb0GF9AsraAL3vIqaTD7vOdma+OVufSxmaom8oyC2tvjlKu5SzBi/80+aUcgPWwdPqv0fx3oyxu/fgfOAFcDFwH+llLZExCBgPpAn8ZsNvLsyPuEoYG1KaWlEXAP8T0S0X2VOBj6U43VURf16HWprg7kfy4qdkGC/M+Hl34MhBXejSSmbGmD6iX2/IjzvtdkYwQVzYd+X9n7/TU/BvB9m8/blSfoAdtkfjnhLNsbt8DfBrj1rlwKy8WgL/gAv/UrPx/VtM5YD4PW/hMvOgYtfDKd+NjtX3Z3np1fDtZ/JYh83A974Wxi/T/44ICuA86ars7/vb9/Kxt3d8IVKUrNLluBseuqZQi1NQ2Gvk+DIt8GeJxRbmGXEznDsf2S3jetg2X2w/L5sbOCmdRBNMGR4lmjtPCM7Bzv0vvW2R5qaYfIR2U0ARMRdPPNx1wSMB8oY3wdwCzAjIqYBjwPnAK8t6bU0QNRGQZaKIqp6th+qkAqj+budFjVKKm9RlSigxa81JQbnzPzaE783XHwzg/vYRbK9kmZzzlhmTh3LzKkDY0qnnrT4jQNekVJ61gCZlFJbRGxzQFBEXEbWcjcuIhaTVeocXNn/u8AcsqkcFpBN5/CmyrpVEfEpsosbwAXthV6kbfrzV7Iv6Ie/EUbumn1J/9U7smIjRX4JX34/rF8OU1/Y92NMPwlGjIfbf9q3xO8f/19WLOTod/U9ho6O/xDc9XOY/Z4s2Rncg+6a98+B6z6TTVEx883FxAFZQvG2a+GXs+DX74Sbv5MlpNNelM05uOXpbML6B+Zk52/LBjhyFpz0MRi6Y3FxQPa+mfGS7PbkP7OEf+kd2VQMrZuyHxXGTIVdD4Y9jnmmYmiZhu0EU47KbqoVHa+HLcATKaWW7jbOI6XUUim8dg1ZknlJSqmYmYNVl4psQy/iWEUUVMlra3GX3IlfMfMs5v2bimzxyz+dQ/55/I6evkcLhp0AACAASURBVDPnHDGZjVtacx1nhyFNHL6HvRF6qidj/LqsVFZZt82B6ymlc7ezPgFdfmtNKV1CVkxGDaL03l1rFsGfPg8HvDyr4hgBzUPhj/+ddcHb57TiXuvhG7P7aTkSv6ZmOOg18PeLsmqLvZnTr3VLluBOORomH9n3GDoaPhbO+iZc8Xr4zb/Dy76z7SqhD10LP39jVqH0jK8W/w88Zo8sAb3jsqw652/f99xtmoZkU2kcc37vWin7asddsykhDsrXiVSNp/OPo/3wenPIfjyVtsoznq2ocXVFdbHMnWy1T+eQM47uCk/0xqCI/MlWW/ux8sXSFFFAV8/8id+4kUP53CsPznUM9V7jVNqQbvxSdv+STz1zxXjBe7L52eZ+PJugu6jJsR/+E4yakrX05PG8c+GmC7NxYUfN6vl+d10Jaxf1fjze9ux3Jpz40WxahrbWLBEc0mlcZkow7xK4+v0wbp9s+oWiu9K2GzQIDn1d1tVzxYOweB5sWAnNw2Dsntk0DGW9tiQNQPkTroJa2fJ2r9za4ldAhdECzkne7pXtuxdRYXTeo6v56twH+3yMpWs3snMdFDLRc5n4qXz9MW5g05PZhOiHnAOjO9Q8aBoML/pP+OXbsnFae5+c/7Xa2uDRv8A+p+c/1q4HZbdbfwBHvLVniWnrlizJnXAAzCjg7+nsRf8Fg5rhD5+EJbfBse/LJhRvGpIlXjd9O/v7p58Er7qkvHFkHUVkY9aKGrsnSQ2k0eolFTqZfO79UwGVTgto8dva1TPXYZixy0hue2wN/3hsTe7jqP6Y+KnfFFlu+jnu/XU2zuvQNzx33f4vy1r8brm4mMTvibuzgiJ5Jxlv94J/zxLT+38D+5+9/e3n/QBWLoBzLy/van/se2H3w7KJ3We/+9nrRu6SdaU97LziWlAlSbnlyS3a9y1m4vRipkAoQjHFXfIpYoxfKmiM3y/feUzOSFTPTPzUGO64HMZOh0ldVBlsHgIHvyar9Ll+ZVYdMY9HKuP78hR26ejAV8KfvgDXfw72eWk29q87Ty2D6z+bvfbepxbz+t3Z8zh4x1+zVr+ld2RdP8fvk40rbNr+fHaSpP7RYA1+hYyreyY/ytnVk/zBBEVM55Dd5x3jp4HNn+tVuiInYu3S02vg0b/CAS/r/ifCg16dzXF276/yv97DN2bjy0ZNzH8syAqonPQxWHZvVr2yOynBb/4DNq+H07/UP317ImDi4VnFziPflrVymvRJUk3Kc7Ut6kodUUwrW+44ikqHi5jOoYBKmmlr4mfmp76zxU/9prSPqoXXZ5NX77WNqal3OTArRHLPVdl8dX3V2pKV8z/wFX0/Rlf2Owv2Pi2bj27qsVmlzM7+8jV44Ldw8qdhwr7Fvr4kSRSTMBUzHUNRxV3yxZFIuVvZBkU2991Dy5/q8zGWPbkRaLzWXfUvEz/VvwVzYdiorrt5touAvU+Bm76TFYLp61xvS26DzU/C9BP6tv+24jvrG/C/J8JPz4HXX5kVfYHsqvXXb2TFVg58JTy/oHn7JEkNodQx9L1UTOKYCuvqmX9sXf4ONjsMbmLjljZO+vKfch0nAoY2b2OaJWk7TPxUuiKrc3Xpoethz+O3PTYOsgm4//oNePiGvk2YDlnrIgFTCyrs0tHICfC6n8OPX54lgAf/C4yaBPN/D4/fms1PePa3LagiSepSnu6EqaCqkXnjyPbPH8PWefwKmEYhbzL7lmP3ZMYuO+aeP2/CjsMYNdzhFuo7Ez/1m1J+kFy7GNYthinv2f62k58PQ0bC/Ln5Er/dDs5fIKY7E/aDWX+Caz+VVSrdtC7ronrmN7KKpSZ9kqQaVti1vrAhevmnUcj7N40aPpgzD9k930GkApj4qb4t+nt2P/nI7W/bPASmHQcL/ti3vhubnspe7+h39j7O3thxFzj7W3Dm17NKms1Dyn09SVJdK6SVLv8hCj1WYV09a2A6B6lW2Hyg0pXa03PxLdC8wzPj4bZnrxNh7WOw8qHev9bDN0DbFtiz4PF93RnUZNInSeqxskdW9EQRSVJKBUya3n6s/OH0TxVtqR+Y+Km+Lbo5m26gp1MMTD8xu194Xe9f64HfwtCdYA8nP5UkqUuRf866gsLIbeu4x/yHkmqCiZ/6TWFz6rTbsjGbWHzyNqp5djZ2Txi9Bzx0be9eq60VHvgdzDjZVjhJUk0pppWtcqwaaN1KFJG4tRd3yVPwpnKk6p8SqRAmfipdab/8Lbs3m5S9qznvtmX6idkk7K1ber7P4ltgwwrY9/TevZYkSf2kJlraKGgKhbxxFDjusfAfrqUqsbiL6tcTd2f3uxzYu/2mnwC3/gAWz4M9ju7ZPvf9BgYN3vYk8ZIkVUMhTVKN1a2x/e+47O+P8beHVvbpGG0FTnEh1QITP/Wbwj84/3l3Nj3DmGm922/aiyCa4MHf9Szxa22Bu36edfMctlPfYpUkqWR5py4oQkQB8/iRv7jLLjsNY7/dduLhFet5eMX6Ph9n91HDOGjSqFyxSLXCxE+lK+1C9MTdMGH/3s9tt8OYbML3e34JL/7k9jPShdfDU0/A887tU5iSJKnniujqOWJoM1ef/8JC4pEahWP8VJ9Sylr8du1lN892B70K1jyWdffcnjt+miWLM07u22tJklSiYou75DuOvSKl2mXip35T6MVgzWOwaW3vx/e12/cMaBoKd/1s29utWwL3/hoOOReah/bttSRJ6g/V7+lZiGKqekrqzMRPpSulytiy+7L7viZ+w3aC/c+G2y+DjWu73+7m70Fqg6P+rW+vI0lSHSiqgmUUMI9ftr+Zn1Q0Ez/Vp5Xzs/txM/p+jKPfCZufhFsv7Xr9k0/ALRdnCeKYqX1/HUmSSlRrrWO1UGRG0nOZ+Kn/FHlhWjEfdhgLw8f2/Ri7H5pV+PzL1+Dp1c9df92noWUjnPixvr+GJEn9pBbSrWIu9anmklmpEZj4qXSlXIhWLsjX2tfu5M9kSd/VH3x2n9T758BtP4Kj3g47T8//OpIk1bDCirtE/iEeRVT1lPRcJn6qTyvmw84FJH67HQwvej/ceTlc8+Gse+cdV8CVb4LdDrG1T5JU8/KOy6s1pdQGkOQ8fuo/hV2YNq6F9ctg3F7FHO+4D8CGlXDTt7MbwMTD4bU/g8HDinkNSZJKlidhah+Xl/dKHeQv7gK1N25RagSlJn4RcSrwdaAJuDil9LlO678KnFB5OhyYkFIaXVnXCtxVWfdYSumsMmNViYr+6W7lguy+iBY/yCaAf+mXsikbFt2UFXLZ+1QY1FTM8SVJKlGjJUmJ1HCtmFItKC3xi4gm4ELgJcBi4JaImJ1Surd9m5TSezts/x7g0A6HeDql9Lyy4lMdW9Ge+BXU4tdu0uHZTZKkOlRENc3cSWTkj8OunlI5yhzjdySwIKW0MKW0GbgcOHsb258LXFZiPKqywn6RXDkfYhCMnVbQASVJGthqLdlqtFZMqRaU2dVzIrCow/PFwFFdbRgRewDTgGs7LB4WEfOAFuBzKaVflRWoylX4tWTlQzB6CjQPLfrIkiTVnVrKkQJYuHw9H77qru1u251bH11dU3+T1CjKTPy6+j/bXQ5wDnBlSqm1w7IpKaUlEbEncG1E3JVSeug5LxIxC5gFMGXKlLwxqx6sXZQlfpIkaatcxV227psv5Tpi6lgeWr6e39/zRK7jvHDGuFz7S3quMhO/xcDkDs8nAUu62fYc4F0dF6SUllTuF0bE9WTj/56T+KWULgIuApg5c2aNdVRQR4X9erd2MUw/qaijSZKkgvznKfvwn6fsU+0wJHWhzDF+twAzImJaRAwhS+5md94oIvYBxgB/67BsTEQMrTweBxwD3Nt5X9WHQscNtGyGJ/8JoyYVeFBJkupXEePhtk7nYB9LqWGV1uKXUmqJiHcD15BN53BJSumeiLgAmJdSak8CzwUuT+lZ6cF+wPcioo0sOf1cx2qgGsDWPQ4kGD15u5tKkjSQ2O1J0raUOo9fSmkOMKfTso93ev7JLvb7K3BQmbGp/0URPyOurdQLssVPkiQA57yT1CNldvWUAEhF9vVcuzi7H2WLnyRJHeW53rbvagopNS4TP9WXNZUWv50mVjcOSZIkqY6Y+KnfFPIr4tpFMGICDB5WxNEkSap/BTbTFTIsQ1JNMvFT6QodbL52kYVdJEnqQqFVtCU1HBM/1Ze1iy3sIkmSJPWSiZ/6Te7eIylVEj9b/CRJaldE50yLu0iNz8RPpSus68n6FdCy0cRPkiRJ6iUTP9WP9jn8HOMnSVKhUmVEvrVdpMZl4qd+k3uCWSdvlyTpOazEKaknTPxUP5y8XZKkbhUxtMIcUmpcJn4qXWHVpdcsgsEjYIcxRR1RkqS6V2RxF0mNy8RP/Sfvlal9Dj9/jpTU4CLiixFxf0TcGRFXRcToDus+FBELIuKBiDilmnGqtqRiZ86V1GBM/FQ/1i5yfJ+kgWIucGBK6WDgQeBDABGxP3AOcABwKvDtiGiqWpRqGO0pY+7x+JJqlomfSpeK6j/i5O2SBoiU0u9TSi2VpzcB7R9+ZwOXp5Q2pZQeBhYAR1YjRtUOO8JI6gkTP/WbXBemzethw0oLu0gaiN4MXF15PBFY1GHd4soyqZhxeiaRUsNqrnYAUo+sfTy7N/GT1CAi4g/Arl2s+khK6deVbT4CtAA/ad+ti+27/LofEbOAWQBTpkzJHa8aW2G9cyTVLBM/1Ye1j2X3Tt4uqUGklF68rfURcR5wBnBSeuZb+WKg4wfhJGBJN8e/CLgIYObMmX6rb2B29ZTUE3b1VL/JdV3aOoefY/wkNb6IOBX4AHBWSmlDh1WzgXMiYmhETANmAH+vRoyqPXmy+2eKu0hqVLb4qT6sXQwxCHbcvdqRSFJ/+BYwFJgbWXPOTSmlt6eU7omInwH3knUBfVdKqbWKcUqS6oSJn0pXyLCBNYuypK/Jt6ykxpdS2msb6z4DfKYfw1GNK3IKhrDfqNSw/Bat+uAcfpIkbdNbLr2FIU19G8Xz9BYbjqVGZ+KnfpPrV8Q1i2DKUcUFI0lSgzhy2ljOPGR3NuVM3g6dPJqDJ44qKCpJtcbET6VLuYabA60tsO5xp3KQJKkLu4/egW+ee2i1w5BU46zqqdr35FJIrTDaeagkSZKkvjDxU7/pc0fPtYuye+fwkyRJkvqk1MQvIk6NiAciYkFEfLCL9W+MiOURcXvl9tYO686LiPmV23llxqly5a7quaaS+I2yxU+SJEnqi9LG+EVEE3Ah8BJgMXBLRMxOKd3badMrUkrv7rTvWOATwEyyOUVvrey7uqx4VcPWPpbdW9VTkiRJ6pMyW/yOBBaklBamlDYDlwNn93DfU4C5KaVVlWRvLnBqSXGqn/S5qOeaRTB8HAwZXmg8kiRJ0kBRZuI3EVjU4fniyrLOXhkRd0bElRHRPoirp/uqDuSev33tIsf3SZIkSTmUmfh11b7TOQf4DTA1pXQw8Afg0l7sm20YMSsi5kXEvOXLl/c5WNWwNY9Z0VOSJEnKoczEbzHQsZlmErCk4wYppZUppU2Vp/8LHN7TfTsc46KU0syU0szx48cXErjKEX2p65kSrF3sHH6SJElSDmUmfrcAMyJiWkQMAc4BZnfcICJ26/D0LOC+yuNrgJMjYkxEjAFOrixTHcpV1XP9cmjZaIufJEmSlENpVT1TSi0R8W6yhK0JuCSldE9EXADMSynNBv49Is4CWoBVwBsr+66KiE+RJY8AF6SUVpUVq2rY1qkcbPGTJEmS+qq0xA8gpTQHmNNp2cc7PP4Q8KFu9r0EuKTM+NS/+lTVs30qB4u7SJIkSX1W6gTuEkDKU9fTFj9JkiQpNxM/1ba1i2DoTrDD6GpHIkmSJNUtEz/VtjWLLOwiSZIk5WTip9Llquq55jG7eUqSJEk5mfipdqUEqx+BMXtUOxJJkiSprpn4qd/0uqrn+uWwZT2MmVZKPJIkSdJAYeKn2rX6kex+rImfJEmSlIeJn2rXqoeze1v8JEmSpFxM/NRvgl729Vz9MBBW9ZQkSZJyMvFT6VJfy3quehh22h0GDys2IEmSJGmAMfFT7Vr9sN08JUmSpAKY+Knf9Lqq5+pHYOzUEiKRJEmSBhYTP5WuTz09N6+Hp56wxU+SJEkqgImfapNTOUiSJEmFMfFTv+lVT8+tUzlMLSESSZIkaWAx8VPp+lTTs73Fz66ekiRJUm4mfqpNqx+GYaNg+NhqRyJJkiTVPRM/9ZvoTVnPVU7lIEmSJBXFxE+1aeV82Hl6taOQJEmSGoKJn0rX6+kctjwNaxbBuL1LiUeSJEkaaEz81G963NFz1UIgwc57lRiNJEmSNHCY+Kn2rJif3Y+bUd04JEmSpAZh4qfSpd5O6LCykvjZ4idJkiQVwsRP/abHRT1XzIedJsKQEaXGI0mSJA0UJn6qPSvm29onSZIkFajUxC8iTo2IByJiQUR8sIv174uIeyPizoj4Y0Ts0WFda0TcXrnNLjNOlatXVT1TgpULrOgpSZIkFai5rANHRBNwIfASYDFwS0TMTind22GzfwAzU0obIuIdwBeAf6msezql9Lyy4lP/69EE7k8tg03rLOwiSZIkFajMFr8jgQUppYUppc3A5cDZHTdIKV2XUtpQeXoTMKnEeFQPLOwiSZIkFa7MxG8isKjD88WVZd15C3B1h+fDImJeRNwUES8rI0D1j17V9FzxYHZvi58kSZJUmNK6etL1fN1d5gAR8XpgJnBch8VTUkpLImJP4NqIuCul9FAX+84CZgFMmTIlf9SqruUPwOARsJONv5IkSVJRymzxWwxM7vB8ErCk80YR8WLgI8BZKaVN7ctTSksq9wuB64FDu3qRlNJFKaWZKaWZ48ePLy56VccT98CEfWGQBWclSZKkopT57foWYEZETIuIIcA5wLOqc0bEocD3yJK+ZR2Wj4mIoZXH44BjgI5FYVRPelrWMyVYdi9M2L/ceCRJkqQBprTEL6XUArwbuAa4D/hZSumeiLggIs6qbPZFYCTw807TNuwHzIuIO4DrgM91qgaqRrR+OWxYCbscUO1IJKkmRMR/RkSq/AhKZL5RmSbpzog4rNoxSpLqQ5lj/EgpzQHmdFr28Q6PX9zNfn8FDiozNvWvnszkwBP3ZPe2+EkSETGZbEqkxzosPg2YUbkdBXynci9J0jY5kEql63FVz2WVRl0TP0kC+Crwfp79MXo28KOUuQkYHRG7VSU6SVJdMfFT7XjiXhgxHkZapEfSwFYZEvF4SumOTqt6O1WSJElAyV09pXY96elpYRdJA0lE/AHYtYtVHwE+DJzc1W5dLOtuqiSnO5IkbWXip9L1qKhnWxssvx8Of2PZ4UhSTehunHtEHARMA+6IbID0JOC2iDiSHk6VVDn+RcBFADNnzuxxr3tJUmOyq6dqw8oFsGUD7HJgtSORpKpKKd2VUpqQUpqaUppKluwdllL6J9m0SP9aqe75fGBtSmlpNeOVJNUHW/zUL2J7ZT2X/CO73/3Q8oORpPo1BzgdWABsAN5U3XAkSfXCxE+lSz2p67n0dhg8HMbtXX5AklRHKq1+7Y8T8K7qRSNJqld29VRtWPIP2PUgaPK3CEmSJKloJn7qF9vs6NnWCkvvsJunJEmSVBITP5Vuu1U9V8zPCrvs9rx+iUeSJEkaaEz8VH1bC7uY+EmSJEllMPFTv9hmUc9FN8PQURZ2kSRJkkpi4qfSbbem52M3weQjYVBTf4QjSZIkDTgmfqquDatg+X0w5fnVjkSSJElqWCZ+6hfRXV3PRTdn91OO7r9gJEmSpAHGxE+l22ZVz8f+BoMGw8TD+i0eSZIkaaAx8VN1LbweJh0Bg3eodiSSJElSwzLxU//oqqfnU8uzidv3OrHfw5EkSZIGEhM/lS51V9fzoWuz++kn9V8wkiRJ0gBk4qfqeeiPMHxn2M2J2yVJkqQymfipXzynp2frFpg/N2vtG+TbUJIkSSqT37hVvq56ei78Ezy9Cg54eb+HI0mSJA00Jn6qjrt/AUNHwV6O75MkSZLKZuKnfhEd+3puXAf3/Qb2OwOah1YtJkmSJGmgMPFT6Z7T0/OOy2HzkzDzLdUIR5IkSRpwSk38IuLUiHggIhZExAe7WD80Iq6orL85IqZ2WPehyvIHIuKUMuNUP9qyEf76TZg4EyYdXu1oJEmSpAGhtMQvIpqAC4HTgP2BcyNi/06bvQVYnVLaC/gq8PnKvvsD5wAHAKcC364cT3Uq2ut63vRtWPsYnPjR6gYkSZIkDSBltvgdCSxIKS1MKW0GLgfO7rTN2cCllcdXAidFRFSWX55S2pRSehhYUDmeqiWlrm9tbV3cWp91i9TKINrg/jlw3WdgvzNh+gnV/oskSZKkAaO5xGNPBBZ1eL4YOKq7bVJKLRGxFti5svymTvtOLC/UzKMXHMD41uU92ja6mKOgq2XZ8p7tv63lZcQwKHr+Wnl8CPhQE1nqv+vBcNY3++V1JUmSJGXKTPy6yjU6ZxrdbdOTfbMDRMwCZlWePhURD/Q4wq6NA1bkPEZ/qcNY/wzvGFvtWHqiDs9tXTDW8tRTvEXEukcRgQwUt95664qIeDTnYerpPQb1Fa+xlqOeYoX6itdYy1FUrF1eI8tM/BYDkzs8nwQs6WabxRHRDIwCVvVwXwBSShcBFxUUMxExL6U0s6jjlclYy1NP8RprOeopVqiveOsp1kaRUhqf9xj19u9WT/EaaznqKVaor3iNtRxlx1rmGL9bgBkRMS0ihpAVa5ndaZvZwHmVx68Crk0ppcrycypVP6cBM4C/lxirJEmSJDWs0lr8KmP23g1cAzQBl6SU7omIC4B5KaXZwPeBH0fEArKWvnMq+94TET8D7gVagHellFrLilWSJEmSGlmZXT1JKc0B5nRa9vEOjzcCr+5m388Anykzvm4U1m20HxhreeopXmMtRz3FCvUVbz3FqmfU279bPcVrrOWop1ihvuI11nKUGmtkPSslSZIkSY2qzDF+kiRJkqQaMGATv4g4NSIeiIgFEfHBLtYPjYgrKutvjoip/R8lRMTkiLguIu6LiHsi4vwutjk+ItZGxO2V28e7OlZ/iIhHIuKuShzzulgfEfGNynm9MyIOq1Kc+3Q4X7dHxLqI+I9O21T1vEbEJRGxLCLu7rBsbETMjYj5lfsx3ex7XmWb+RFxXlfb9EOsX4yI+yv/zldFxOhu9t3me6afYv1kRDze4d/69G723ebnRj/Ge0WHWB+JiNu72be/z22Xn1e1+r5V17w+lsdrZGHx1c31cRvxeo0sJ1avj9uSUhpwN7JiMw8BewJDgDuA/Ttt807gu5XH5wBXVCnW3YDDKo93BB7sItbjgf+r9nmtxPIIMG4b608Hriabq/H5wM01EHMT8E9gj1o6r8CLgMOAuzss+wLwwcrjDwKf72K/scDCyv2YyuMxVYj1ZKC58vjzXcXak/dMP8X6SeA/e/A+2ebnRn/F22n9l4GP18i57fLzqlbft966/Df0+lhuzF4ji4mpbq6P24jXa2QJsXZa7/Wx022gtvgdCSxIKS1MKW0GLgfO7rTN2cCllcdXAidFRFcTy5cqpbQ0pXRb5fGTwH3AxP6Oo0BnAz9KmZuA0RGxW5VjOgl4KKWUd3LjQqWUbiCrdttRx/flpcDLutj1FGBuSmlVSmk1MBc4tbRA6TrWlNLvU0otlac3kc3HWXXdnNee6MnnRuG2FW/lM+k1wGVlx9ET2/i8qsn3rbrk9bG6vEb2QD1dH8FrZFm8PvbeQE38JgKLOjxfzHMvFlu3qfzHXAvs3C/RdaPSneZQ4OYuVh8dEXdExNURcUC/BvZsCfh9RNwaEbO6WN+Tc9/fzqH7D4ZaOa/tdkkpLYXsQwSY0MU2tXiO30z2K3ZXtvee6S/vrnS5uaSbrha1eF5fCDyRUprfzfqqndtOn1f1+r4diLw+lstrZHnq+XPGa2TxvD52YaAmfl39Mtm5vGlPtuk3ETES+AXwHymldZ1W30bWBeMQ4JvAr/o7vg6OSSkdBpwGvCsiXtRpfa2d1yHAWcDPu1hdS+e1N2rtHH+EbD7On3SzyfbeM/3hO8B04HnAUrLuIZ3V1HmtOJdt/5pZlXO7nc+rbnfrYlm1z+9A5PWxXF4jq6umzi94jSyR18cuDNTEbzEwucPzScCS7raJiGZgFH1r+s4tIgaTvUl+klL6Zef1KaV1KaWnKo/nAIMjYlw/h9key5LK/TLgKrKm/456cu7702nAbSmlJzqvqKXz2sET7d1+KvfLutimZs5xZQDyGcDrUqWjemc9eM+ULqX0REqpNaXUBvxvNzHUzHmFrZ9LrwCu6G6bapzbbj6v6up9O8B5fSyR18hS1d3njNfIcnh97N5ATfxuAWZExLTKr1nnALM7bTMbaK+a8yrg2u7+U5ap0kf5+8B9KaWvdLPNru3jKyLiSLJ/15X9F+XWOEZExI7tj8kGLt/dabPZwL9G5vnA2vYm7irp9hehWjmvnXR8X54H/LqLba4BTo6IMZXuGCdXlvWriDgV+ABwVkppQzfb9OQ9U7pOY2he3k0MPfnc6E8vBu5PKS3uamU1zu02Pq/q5n0rr49l8RpZurr6nPEaWSqvj91J/VTNptZuZJWzHiSrQPSRyrILyP4DAgwj69qwAPg7sGeV4jyWrDn3TuD2yu104O3A2yvbvBu4h6yC0k3AC6oU656VGO6oxNN+XjvGGsCFlfN+FzCziu+B4WQXqVEdltXMeSW72C4FtpD92vMWsnE0fwTmV+7HVradCVzcYd83V967C4A3VSnWBWR90tvft+1VAHcH5mzrPVOFWH9ceT/eSfYhvFvnWCvPn/O5UY14K8t/2P5e7bBttc9td59XNfm+9dbtv6PXx3Li9RpZXGx1c33cRrxeI0uItbL8h3h97PIWlYNJkiRJkhrUQO3qKUmSJEkDhomf8g1UDgAAAD1JREFUJEmSJDU4Ez9JkiRJanAmfpIkSZLU4Ez8JEmSJKnBmfhJkiRJUoMz8ZMkSZKkBmfiJ0mSJEkN7v8HKaFGVCagx9QAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1080x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_outputs(result)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}